if (!require('knitr')) install.packages('knitr'); library('knitr')
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.align='center')
# Load in packages
if (!require("pacman")) install.packages("pacman"); library(pacman) # for rapid install if not in library
# load packages
pacman::p_load("RgoogleMaps", "SDMTools", "rgdal", "ggmap", "gridExtra", "automap", "cowplot","sp", "gstat", "plotrix", "dplyr", "ggplot2", "scales","magrittr","plyr","effects", "ggpubr", "sjstats", "RColorBrewer")
Hakalau National Forest Wildlife Refuge Isoscapes
Data collected 20-21 August 2019 in Hakalau Forest, Hawai’i Island. Soil and plant samples collected from afforested/restored plots after a century of deforestation from cattle grazing (AK sites) and remnant mixed koa-ohia forests (RK sites). Target engineers of forest habitats are Acacia koa (Hawaiian name: Koa), an indigenous hardwood species that forms native forest canopies and associates with nitrogen fixing bacteria. Understory species include Rubus argutus – an invasive species (common name: sawtooth blackberry) – and an endemic species Rubus hawaiiensis (common name: ’Ākala); both are members of the Roseaceae family. Noteably, R. argutus was found more in the remnant plot (i.e., RK) and the indigenous R. hawaiiensis was common in the restored plot (i.e, AK).
Soil samples were collected in a spatially explicit fashion (every 4 or 5m) with a 20 x 35m superplot, consisting of thirty-five 4 x 5m subplots. This data was used to create a δ13C and δ15N stable isotope landscape, hereafter ‘isoscape’. Koa and Rubus spp. were collected within each 4 x 5m sublot.
Plot densities were…
AK: Koa (18 trees) Rubus (14 plants) per 20 x 35m plot= 0.026 Koa m-2, 0.020 Rubus m-2 RK: Koa (10 trees) Rubus (28 plants) per 20 x 35m plot= 0.014 Koa m-2, 0.040 Rubus m-2
# load data
Hak.gps<-read.csv("data/Hakalau.gps.csv")
API<-read.csv("data/API_key.csv")
API.key<-API[1,1]
ggplot(Hak.gps, aes(x = longitude, y = latitude)) +
coord_quickmap() +
geom_point()
Figure 1. Site map of (a) the Hawaiian Island archipelago, (b) sampling plot in Hakalau Forest Refuge, and (c) plot layouts, highlighting Acacia koa density.
######## using ggmap
register_google(key=API.key)
########
#Hawaiian islands Map
########
hi_map<-get_map(location=c(-157,20.5),zoom=7,maptype="satellite",color="color")
hi_map_for_man <- ggmap(hi_map) +
geom_point(aes(x = -155.320, y = 19.83), pch=23,colour="black",fill="red", size = 2, stroke=0.5) +
xlab("") + ylab("Latitude") +
scale_y_continuous(limits=c(18.8, 22.2))+
scale_x_continuous(limits=c(-160.2, -154)) +
theme(axis.text=element_text(colour="black",size=5),
axis.title=element_text(colour="black",size=8),
plot.margin = unit(c(0, 0.5, 0, 0.7), "cm")) +
ggsn::scalebar(x.min=-160.2, x.max=-154, y.min=18.8, y.max=22.2, dist=50, dist_unit="km", transform=TRUE,
st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84",st.color="white", border.size=0.5)
##########
# site map plot showing distance from plots
##########
Hak.rev=c(-155.317, 19.82158002)
map2<-get_map(Hak.rev,
zoom=16,
scale = 2,
maptype= "satellite",
source="google", extent= "device", legend="topright")
## site single point only
site.only<-Hak.gps[c(1,5),]
Hakalau.site<-
ggmap(map2)+
geom_point(aes(x=longitude, y=latitude), data=site.only, alpha=0.5, color="white", size=4)+
labs(x="", y="") +
scale_y_continuous(limits=c(19.81785, 19.8242))+
scale_x_continuous(limits=c(-155.3223, -155.311)) +
theme(text = element_text(size=6),
plot.margin = unit(c(0.2, 0.5, 0.2, 0.2), "cm")) +
annotate("text", x=-155.3187, y=19.8225, label= "AK", size=3, col="white") +
annotate("text", x=-155.31478, y=19.8224, label= "RK", size=3, col="white") +
ggsn::scalebar(x.min=-155.314, x.max=-155.311, y.min=19.818, y.max=19.824, dist=100,
dist_unit="m", transform=TRUE,
st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84",st.color="white", border.size=0.5)
Hakalau.site
Figure 1. Site map of (a) the Hawaiian Island archipelago, (b) sampling plot in Hakalau Forest Refuge, and (c) plot layouts, highlighting Acacia koa density.
##########
# Site plots: AK or RK only
##########
AK<- Hak.gps[(Hak.gps$Site=="AK"),]
RK<- Hak.gps[(Hak.gps$Site=="RK"),]
####### AK map
AK.gps<-c(-155.3189, lat=19.8219)
mapAK<-get_map(AK.gps,
zoom=19,
scale = 2,
maptype= "satellite",
source="google", extent= "device", legend="topright")
AK.map<-ggmap(mapAK, group=type)+
geom_point(aes(x=longitude, y=latitude, shape=type, color=type), data=AK, alpha=0.5, size=3)+
scale_y_continuous(limits=c(19.82165, 19.82225))+
scale_x_continuous(limits=c(-155.3194, -155.3183)) +
theme(legend.position="top",
legend.text = element_text(color = "white"),
legend.title = element_text(color = "white"),
legend.key = element_rect(fill = "white"),
plot.margin = unit(c(0, 0.5, 0, 0), "cm"),
text = element_text(size=8)) +
guides(colour= guide_legend(override.aes = list(color = "white"))) +
scale_color_manual(values=c('red','mediumseagreen'))+
labs(x="Longitude", y="Latitude") +
annotate("text", x=-155.31925, y=19.82225, label= "Restored Forest (AK)", size=3, col="white")+
ggsn::scalebar(x.min=-155.3183, x.max=-155.3187, y.min=19.82165, y.max=19.82225, dist=10,
dist_unit="m", transform=TRUE,
st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84", st.color="white", border.size=0.5)
####### RK map
RK.edit<-RK[-12,] # remove "proximate" Koa
RK.edit[9,3]<-"koa" # rename to be "koa"
RK.gps<-c(-155.315, 19.8216)
mapRK<-get_map(RK.gps,
zoom=19,
scale = 2,
maptype= "satellite",
source="google", extent= "device", legend="topright")
RK.map<-ggmap(mapRK, group=type)+
geom_point(aes(x=longitude, y=latitude, shape=type, color=type), data=RK.edit, alpha=0.5, size=3)+
scale_color_manual(values=c('red','mediumseagreen', "dodgerblue"))+
scale_shape_manual(values=c(16, 17, 3)) +
scale_y_continuous(limits=c(19.82120, 19.82195))+
scale_x_continuous(limits=c(-155.31575, -155.3144)) +
theme(legend.position="top",
legend.title=element_blank(),
legend.key = element_rect(fill = "white"),
plot.margin = unit(c(0, 0.5, 0, 0), "cm"),
text = element_text(size=8)) +
labs(x="Longitude", y="") +
annotate("text", x=-155.3156, y=19.82193, label= "Remnant Forest (RK)", size=3, col="white")+
ggsn::scalebar(x.min=-155.3148, x.max=-155.3144, y.min=19.8212, y.max=19.8219, dist=10,
dist_unit="m", transform=TRUE,
st.bottom=FALSE, st.size=2, box.fill=c("gray50", "white"), model="WGS84", st.color="white", border.size=0.5)
site.plots<-plot_grid(hi_map_for_man, Hakalau.site, AK.map, RK.map,
labels=c('a', 'b', 'c', 'd'), label_size=8, hjust=-1, vjust= 2, ncol=2, nrow=2)
### export it
pdf(file= "figures/Fig1.sitemaps.pdf", height=7, width=8)
site.plots
dev.off()
## quartz_off_screen
## 2
Hak.df<-read.csv("data/Hakalau.isotopes.csv")
########## DBH
DBH.df<- Hak.df %>%
select(c(DBH..cm.1, DBH..cm.2, DBH..cm.3, DBH..cm.4, DBH..cm.5))
# replace NAs with zeros
DBH.df<- DBH.df %>%
mutate_each(funs(replace(., which(is.na(.)), 0)))
# if only one trunk, then make new column and report this value
DBH.df <- DBH.df %>%
mutate(DBH.1trunk = ifelse(DBH..cm.2>"0", 0, DBH..cm.1))
# if multiple trunks, sum of squares and square root all trunks
# and plug in the DBH for single trunks
DBH.df <- DBH.df %>%
mutate(DBH.Total = ifelse(DBH.1trunk=="0",
sqrt(DBH..cm.1^2 + DBH..cm.2^2 + DBH..cm.3^2 + DBH..cm.4^2 + DBH..cm.5^2), DBH.1trunk))
# replace zeros with NA now that total DBH has been calculated
DBH.df$DBH.Total<-as.numeric(ifelse(DBH.df$DBH.Total=="0", "NA", DBH.df$DBH.Total))
############ Canopy diameter
# use area of an elipse, a x b x π, where a and b are major and minor radius
Can.df<-Hak.df %>%
select(c(canopy.0..m, canopy.90..m, canopy.180..m, canopy.270..m))
# Major radius is the average of 0 and 180 degrees
Can.df<- Can.df %>%
mutate(Maj.rad = rowMeans(cbind(canopy.0..m, canopy.180..m), na.rm=T))
# Minor radius is the average of 90 and 270 degrees
Can.df<- Can.df %>%
mutate(Min.rad = rowMeans(cbind(canopy.90..m, canopy.270..m), na.rm=T))
# Elipse area
Can.df<- Can.df %>%
mutate(Canopy.area = ifelse(Maj.rad>0, Maj.rad*Min.rad*3.14, "NA"))
########################
# combine DBH.Total and Canopy.area with isotope dataframe and reduce columns
Hak.df$DBH.Total<-DBH.df$DBH.Total
Hak.df$Canopy.area<-Can.df$Canopy.area
# Hak.data is the full dataframe
Hak.data<- Hak.df %>%
select(c(Plot, Position.point, Sample, DBH.Total, Canopy.area, d15N, d13C, N..percent, Total.N..mmol.gdw, C..percent, Total.C..mmol.gdw, C.N))
########################
########################
########################
# test plots
### Scatter of d13C vs. d15N
iso.scatter<-ggplot(data=Hak.data, aes(x=d13C, y=d15N, color=Plot, shape=Sample))+
geom_point(cex=1.5) +
xlab(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))) +
ylab(expression(paste(delta^{15}, N, " (\u2030, air)")))
### Scatter of total C vs. total N
total.nut<-ggplot(data=Hak.data, aes(x=Total.C..mmol.gdw, y=Total.N..mmol.gdw, color=Plot, shape=Sample))+
geom_point(cex=1.5) +
xlab("total Carbon") +
ylab("total Nitrogen")
# figures together
plot_grid(iso.scatter, total.nut,
labels=c('a', 'b'), label_size=10, hjust=-1, vjust= 3, ncol=2, nrow=1)
Ndfa<-read.csv("data/Ndfa.Rub.csv")
# rename the d15N so we know it is for rubus
Ndfa$d15N.Rub<-Ndfa$d15N
Ndfa$d15N.Koa<-Ndfa$mean.d15N.Ndfa
####################
#################### relationship between koa and rubus d15N?
# Yes, but a bit stronger in the RK site
mod.AK<-lm(Ndfa[(Ndfa$Plot=="AK"),]$d15N.Rub~Ndfa[(Ndfa$Plot=="AK"),]$d15N.Koa); anova(mod.AK)
## Analysis of Variance Table
##
## Response: Ndfa[(Ndfa$Plot == "AK"), ]$d15N.Rub
## Df Sum Sq Mean Sq F value Pr(>F)
## Ndfa[(Ndfa$Plot == "AK"), ]$d15N.Koa 1 1.9786 1.97856 4.4696 0.05612 .
## Residuals 12 5.3120 0.44267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
int.rubAK<-coef(summary(mod.AK))[1] #intercept
slop.rubAK<-coef(summary(mod.AK))[2] #slope
mod.RK<-lm(Ndfa[(Ndfa$Plot=="RK"),]$d15N.Rub~Ndfa[(Ndfa$Plot=="RK"),]$d15N.Koa); anova(mod.RK)
## Analysis of Variance Table
##
## Response: Ndfa[(Ndfa$Plot == "RK"), ]$d15N.Rub
## Df Sum Sq Mean Sq F value Pr(>F)
## Ndfa[(Ndfa$Plot == "RK"), ]$d15N.Koa 1 5.503 5.5030 8.9346 0.006043 **
## Residuals 26 16.014 0.6159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
int.RKoa<-coef(summary(mod.RK))[1] #intercept
slop.RKoa<-coef(summary(mod.RK))[2] #slope
#### plot
BW.back<-theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black", size=0.2))
Rub.Koa.d15N.plot<-ggplot(data=Ndfa, aes(x=d15N.Koa, y=d15N.Rub, color=Plot))+
geom_point(aes(color=Plot), alpha=.8) +
scale_color_manual(values=c('gray30', "gray70")) +
coord_equal() + theme_bw() +
xlab(expression(paste(italic("Acacia koa "), delta^{15}, N, " (\u2030)"), sep=""))+
ylab(expression(paste(italic("Rubus"), " spp. ", delta^{15}, N, " (\u2030)"), sep="")) +
guides(colour=guide_legend(override.aes = list(stroke=1, fill=NA, alpha=1))) +
geom_smooth(method = "lm", alpha = .15, size=0.6) +
theme(legend.title = element_blank(),
axis.text=element_text(size=6),
axis.title=element_text(size=10))+
annotate(geom="text", x=2.7, y=0.8, label=expression(paste("AK ", italic("p"), "=0.056")), size=2) +
annotate(geom="text", x=2.7, y=0.6, label=expression(paste("RK ", italic("p="), bold("0.006"))), size=2) +
BW.back
Rub.Koa.d15N.plot
Figure 2a. Relationship between d15N-Koa and 15N-Rubus at two sampling sites.
dev.copy(pdf, "figures/Fig 3.Koa.Rub.regr.pdf", width=4, height=4, encod="MacRoman")
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
# ** No relationship between d15N and DBH or Canopy at plot level**
# ** Effects at PLOTs on d15N (lower in AK) **
######## plot of d15N and DBH, and d15N and Canopy area
# all samples
par(mfrow=c(1,2), mar=c(4,5,1,1))
#### DBH
mod.DBH<-lm(d15N~DBH.Total*Plot, data=Hak.data) # no interaction effect
mod.DBH<-lm(d15N~DBH.Total+Plot, data=Hak.data)
int<-coef(summary(mod.DBH))[1]
DBH.slope<-coef(summary(mod.DBH))["DBH.Total",1] # DBH (slope for all figures)
RK.int<-coef(summary(mod.DBH))["PlotRK",1] # RK, AK = 0
# no significant relationship to d15N but a trend for DECREASE with DBH
plot(d15N~DBH.Total, data=Hak.data,
col=c("coral","dodgerblue")[as.factor(Plot)],
pch=c(1,2)[as.factor(Plot)],
ylab=expression(paste(delta^{15}, N, " (\u2030, air)")),
xlab="DBH (cm)",
ylim=c(0,5), xlim=c(0,60))
ablineclip(int, DBH.slope, col="coral", lwd=2,
x1 = min(Hak.data$DBH.Total[Hak.data$Plot=="AK"], na.rm=T),
x2 = max(Hak.data$DBH.Total[Hak.data$Plot=="AK"], na.rm=T)) # AK model
ablineclip((int+ RK.int), DBH.slope, col="dodgerblue", lwd=2,
x1 = min(Hak.data$DBH.Total[Hak.data$Plot=="RK"], na.rm=T),
x2 = max(Hak.data$DBH.Total[Hak.data$Plot=="RK"], na.rm=T)) # RK model
#### Canopy
mod.Can<-lm(d15N~Canopy.area*Plot, data=Hak.data) # no interaction
mod.Can<-lm(d15N~Canopy.area+Plot, data=Hak.data)
int<-coef(summary(mod.Can))[1]
Can.slope<-coef(summary(mod.Can))["Canopy.area",1] # DBH (slope for all figures)
RK.int<-coef(summary(mod.Can))["PlotRK",1] # RK, AK = 0
# no significant relationship to d15N but a trend for INCREASE with canopy area
plot(d15N~Canopy.area, data=Hak.data, yaxt="n",
col=c("coral","dodgerblue")[as.factor(Plot)],
pch=c(1,2)[as.factor(Plot)],
xlab=expression(paste("Canopy (m"^2,")")),
ylab="",
ylim=c(0,5), xlim=c(0,80))
axis(side=2, labels=F)
ablineclip(int, Can.slope, col="coral", lwd=2,
x1 = min(Hak.data$Canopy.area[Hak.data$Plot=="AK"], na.rm=T),
x2 = max(Hak.data$Canopy.area[Hak.data$Plot=="AK"], na.rm=T)) # AK model
ablineclip((int+ RK.int), Can.slope, col="dodgerblue", lwd=2,
x1 = min(Hak.data$Canopy.area[Hak.data$Plot=="RK"], na.rm=T),
x2 = max(Hak.data$Canopy.area[Hak.data$Plot=="RK"], na.rm=T)) # RK model
legend("topright", legend=levels(Hak.data$Plot), box.lty=0, bg="transparent", lty=1, pch=c(1,2),
col=c("coral", "dodgerblue"), cex=1)
dev.copy(pdf, "figures/d15N.DBH.Canopy.pdf", width=7, height=5, encod="MacRoman")
dev.off()
With Plot being most important effect, divide DBH and Canopy to be plot mean +/- SE
Supplemental Figure 1. with afforested and remnant forests
- DBH and Canopy figure as bar chart and means - using Mann-Whitney nonparametric tests: significant plot effects - greater DBH (AK), trend for greater Canopy (RK)
## Models
mod.DBH<-lm(DBH.Total~Plot, data=Hak.data); anova(mod.DBH) # signif
## Analysis of Variance Table
##
## Response: DBH.Total
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 1290.6 1290.58 8.369 0.00762 **
## Residuals 26 4009.5 154.21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.Can<-lm(Canopy.area~Plot, data=Hak.data); anova(mod.Can) # not signif, a bit wobbly in assumptions
## Analysis of Variance Table
##
## Response: Canopy.area
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 654.2 654.16 2.0862 0.1606
## Residuals 26 8152.6 313.56
# diagnostic plots
for(i in c(4:5)){
Y<-Hak.data[,i]
full<-lm(Y~Plot, data=Hak.data, na.action=na.exclude)
R <- resid(full) #save glm residuals
op<-par(mfrow = c(2,3), mar=c(5,4,1,2), pty="sq")
plot(full, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = colnames(Hak.data)[i])
QQline <- qqline(R)
hist(R, xlab="Residuals", main = colnames(Hak.data)[i])
plot(Hak.data$Plot, R, xlab=colnames(Hak.data)[i], ylab="Residuals")
}
# Mann-Whitney U-Test
mwu(Hak.data, DBH.Total, Plot) # signif
##
## # Mann-Whitney-U-Test
##
## Groups 1 = AK (n = 18) | 2 = RK (n = 10):
## U = 315.000, W = 144.000, p = 0.010, Z = 2.589
## effect-size r = 0.489
## rank-mean(1) = 17.50
## rank-mean(2) = 9.10
mwu(Hak.data, Canopy.area, Plot) # not
##
## # Mann-Whitney-U-Test
##
## Groups 1 = AK (n = 18) | 2 = RK (n = 10):
## U = 252.000, W = 81.000, p = 0.666, Z = -0.432
## effect-size r = 0.082
## rank-mean(1) = 14.00
## rank-mean(2) = 15.40
#### plots
DBH.m<-aggregate(DBH.Total~Plot, data=Hak.data, mean)
Can.m<-aggregate(Canopy.area~Plot, data=Hak.data, mean)
DBH.SE<-aggregate(DBH.Total~Plot, data=Hak.data, std.error)
Can.SE<-aggregate(Canopy.area~Plot, data=Hak.data, std.error)
DBH.n<-aggregate(DBH.Total~Plot, data=Hak.data, length)
Can.n<-aggregate(Canopy.area~Plot, data=Hak.data, length)
mean.eco<-cbind(DBH.m, Can.m[2], DBH.SE[2], Can.SE[2]); colnames(mean.eco)<-c("Plot", "DBH", "Can", "D.SE", "C.SE")
pd <- position_dodge(0.7) #offset for error bars
BW.back2<-theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black", size=0.4))
### DBH plot of mean/SE by plot
DBH.plot<-ggplot(data=mean.eco, aes(x=Plot, y=DBH))+
geom_bar(stat="identity", position = pd, width=0.7, fill="gray70")+
geom_errorbar(aes(ymin=DBH-D.SE, ymax=DBH+D.SE), size=.5, width=0, position=pd) +
xlab("Plots") +
ylab("dbh (cm)") +
theme(legend.title = element_blank(),
axis.text=element_text(size=6),
axis.title=element_text(size=8))+
annotate(geom="text", x=1.5, y=35, label=expression(paste(bold("*"))), size=4) +
BW.back2
### Canopy plot of mean/SE by plot
Canopy.plot<-ggplot(data=mean.eco, aes(x=Plot, y=Can))+
geom_bar(stat="identity", position = pd, width=0.7, fill="gray70")+
geom_errorbar(aes(ymin=Can-C.SE, ymax=Can+C.SE), size=.5, width=0, position=pd) +
xlab("Plots") +
ylab(expression(paste("Canopy (m"^2,")"))) +
theme(legend.title = element_blank(),
axis.text=element_text(size=6),
axis.title=element_text(size=8))+
BW.back2
# figures together
plot_grid(DBH.plot, Canopy.plot,
labels=c('a', 'b'), label_size=8, hjust=-1, vjust= 3, ncol=2, nrow=1)
dev.copy(pdf, "figures/Supp.Fig1.DBHCanopy.pdf", width=4, height=3.5, encod="MacRoman")
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
Run statistics on the responses (dd15N, 13C, tissue C and N)
# models
## tests for assumptions
for(i in c(6:12)){
Y<-Hak.data[,i]
full<-lm(Y~Plot*Sample, data=Hak.data, na.action=na.exclude)
R <- resid(full) #save glm residuals
op<-par(mfrow = c(2,3), mar=c(5,4,1,2), pty="sq")
plot(full, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = colnames(Hak.data)[i])
QQline <- qqline(R)
hist(R, xlab="Residuals", main = colnames(Hak.data)[i])
plot(Hak.data$Plot, R, xlab=colnames(Hak.data)[i], ylab="Residuals")
plot(Hak.data$Sample, R, xlab=colnames(Hak.data)[i], ylab="Residuals")
}
########## Separate dataframes
soil<-Hak.data[(Hak.data$Sample=="Soil"),]
Koa<-Hak.data[(Hak.data$Sample=="Koa"),]
Rubus<-Hak.data[(Hak.data$Sample=="RUBARG" | Hak.data$Sample=="RUBHAW"),]
AK.df<-Hak.data[(Hak.data$Plot=="AK"),]; AK.df.sansoil<-AK.df[!(AK.df$Sample=="Soil"),]
RK.df<-Hak.data[(Hak.data$Plot=="RK"),]; RK.df.sansoil<-RK.df[!(RK.df$Sample=="Soil"),]
d15N model
######## Nitrogen isotopes
d15N.soil<-lm(d15N~Plot, data=soil); anova(d15N.soil) # soil signif.
## Analysis of Variance Table
##
## Response: d15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 2.996 2.9963 2.479 0.1187
## Residuals 94 113.612 1.2086
plot(allEffects(d15N.soil))
d15N.Koa<-lm(d15N~Plot, data=Koa); anova(d15N.Koa) # Koa signif.
## Analysis of Variance Table
##
## Response: d15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 4.7952 4.7952 7.4173 0.01139 *
## Residuals 26 16.8088 0.6465
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d15N.Koa))
d15N.Rubus<-lm(d15N~Plot, data=Rubus); anova(d15N.Rubus) # Rubus signif.
## Analysis of Variance Table
##
## Response: d15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 9.9155 9.9155 13.768 0.0006289 ***
## Residuals 40 28.8076 0.7202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d15N.Rubus))
d15N.AK<-lm(d15N~Sample, data=AK.df.sansoil); anova(d15N.AK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: d15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 2.9463 2.94632 6.0621 0.01978 *
## Residuals 30 14.5807 0.48602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d15N.AK))
d15N.RK<-lm(d15N~Sample, data=RK.df.sansoil); anova(d15N.RK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: d15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 4.4682 4.4682 5.1829 0.02886 *
## Residuals 36 31.0357 0.8621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d15N.RK))
d13C model
######## ########
######## Carbon isotopes
d13C.soil<-lm(d13C~Plot, data=soil); anova(d13C.soil) # soil signif.
## Analysis of Variance Table
##
## Response: d13C
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 24.13 24.1302 22.544 7.342e-06 ***
## Residuals 94 100.62 1.0704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d13C.soil))
d13C.Koa<-lm(d13C~Plot, data=Koa); anova(d13C.Koa) # Koa signif.
## Analysis of Variance Table
##
## Response: d13C
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 12.825 12.8250 21.018 0.0001006 ***
## Residuals 26 15.865 0.6102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d13C.Koa))
d13C.Rubus<-lm(d13C~Plot, data=Rubus); anova(d13C.Rubus) # Rubus signif.
## Analysis of Variance Table
##
## Response: d13C
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 21.848 21.8484 46.101 3.684e-08 ***
## Residuals 40 18.957 0.4739
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d13C.Rubus))
d13C.AK<-lm(d13C~Sample, data=AK.df.sansoil); anova(d13C.AK) #only koa vs. rubus comparision
## Analysis of Variance Table
##
## Response: d13C
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 0.851 0.85100 1.2075 0.2806
## Residuals 30 21.142 0.70474
plot(allEffects(d13C.AK))
d13C.RK<-lm(d13C~Sample, data=RK.df.sansoil); anova(d13C.RK) #only koa vs. rubus comparision
## Analysis of Variance Table
##
## Response: d13C
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 1.4676 1.46758 3.8622 0.05714 .
## Residuals 36 13.6796 0.37999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(d13C.RK))
Nitrogen content model
####### Nitrogen content
TN.soil<-lm(Total.N..mmol.gdw~Plot, data=soil); anova(TN.soil) # soil signif.
## Analysis of Variance Table
##
## Response: Total.N..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.1087 0.108676 1.6099 0.2076
## Residuals 94 6.3456 0.067506
plot(allEffects(TN.soil))
TN.Koa<-lm(Total.N..mmol.gdw~Plot, data=Koa); anova(TN.Koa) # Koa signif.
## Analysis of Variance Table
##
## Response: Total.N..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.03271 0.032711 0.5942 0.4478
## Residuals 26 1.43136 0.055052
plot(allEffects(TN.Koa))
TN.Rubus<-lm(Total.N..mmol.gdw~Plot, data=Rubus); anova(TN.Rubus) # Rubus signif.
## Analysis of Variance Table
##
## Response: Total.N..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.2110 0.21100 3.6973 0.06164 .
## Residuals 40 2.2828 0.05707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(TN.Rubus))
N.AK<-lm(Total.N..mmol.gdw~Sample, data=AK.df.sansoil); anova(N.AK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: Total.N..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 1.5930 1.5930 36.455 1.254e-06 ***
## Residuals 30 1.3109 0.0437
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(N.AK))
N.RK<-lm(Total.N..mmol.gdw~Sample, data=RK.df.sansoil); anova(N.RK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: Total.N..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 0.38328 0.38328 5.7415 0.02189 *
## Residuals 36 2.40323 0.06676
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(N.RK))
Carbon content model
######## ########
######## Carbon content
TC.soil<-lm(Total.C..mmol.gdw~Plot, data=soil); anova(TC.soil) # soil signif.
## Analysis of Variance Table
##
## Response: Total.C..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 81.68 81.678 4.9447 0.02857 *
## Residuals 94 1552.71 16.518
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(TC.soil))
TC.Koa<-lm(Total.C..mmol.gdw~Plot, data=Koa); anova(TC.Koa) # Koa signif.
## Analysis of Variance Table
##
## Response: Total.C..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 0.0377 0.03768 0.0962 0.7589
## Residuals 26 10.1849 0.39173
plot(allEffects(TC.Koa))
TC.Rubus<-lm(Total.C..mmol.gdw~Plot, data=Rubus); anova(TC.Rubus) # Rubus signif.
## Analysis of Variance Table
##
## Response: Total.C..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 34.074 34.074 8.6486 0.005419 **
## Residuals 40 157.595 3.940
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(TC.Rubus))
C.AK<-lm(Total.C..mmol.gdw~Sample, data=AK.df.sansoil); anova(C.AK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: Total.C..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 122.081 122.081 473.95 < 2.2e-16 ***
## Residuals 30 7.728 0.258
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(C.AK))
C.RK<-lm(Total.C..mmol.gdw~Sample, data=RK.df.sansoil); anova(C.RK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: Total.C..mmol.gdw
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 32.592 32.592 7.3308 0.0103 *
## Residuals 36 160.052 4.446
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(C.RK))
C:N
####### C:Ncontent
CN.soil<-lm(C.N~Plot, data=soil); anova(CN.soil) # not signif.
## Analysis of Variance Table
##
## Response: C.N
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 3.349 3.3488 1.7346 0.191
## Residuals 94 181.479 1.9306
CN.Koa<-lm(C.N~Plot, data=Koa); anova(CN.Koa) # not signif.
## Analysis of Variance Table
##
## Response: C.N
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 2.289 2.2886 0.3456 0.5617
## Residuals 26 172.171 6.6220
CN.Rubus<-lm(C.N~Plot, data=Rubus); anova(CN.Rubus) # not signif
## Analysis of Variance Table
##
## Response: C.N
## Df Sum Sq Mean Sq F value Pr(>F)
## Plot 1 4.27 4.2660 0.5037 0.482
## Residuals 40 338.78 8.4694
CN.AK<-lm(C.N~Sample, data=AK.df.sansoil); anova(CN.AK) #only koa vs. rubus comparision, signif.
## Analysis of Variance Table
##
## Response: C.N
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 75.46 75.460 12.771 0.001214 **
## Residuals 30 177.26 5.909
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(allEffects(CN.AK))
CN.RK<-lm(C.N~Sample, data=RK.df.sansoil); anova(CN.RK) #only koa vs. rubus comparision, not signif
## Analysis of Variance Table
##
## Response: C.N
## Df Sum Sq Mean Sq F value Pr(>F)
## Sample 1 24.48 24.482 2.6413 0.1128
## Residuals 36 333.68 9.269
plot(allEffects(CN.RK))
Figure 2
data means and SE for plots
- bar plots of d13C, d15N, TC and TN for all 4 samples types with asterisks for significant effects
## make some means and SE
d13C.m<-aggregate(d13C~Plot+Sample, data=Hak.data, mean)
d15N.m<-aggregate(d15N~Plot+Sample, data=Hak.data, mean)
TN.m<-aggregate(Total.N..mmol.gdw~Plot+Sample, data=Hak.data, mean)
TC.m<-aggregate(Total.C..mmol.gdw~Plot+Sample, data=Hak.data, mean)
CN.m<-aggregate(C.N~Plot+Sample, data=Hak.data, mean)
d13C.SE<-aggregate(d13C~Plot+Sample, data=Hak.data, std.error)
d15N.SE<-aggregate(d15N~Plot+Sample, data=Hak.data, std.error)
TN.SE<-aggregate(Total.N..mmol.gdw~Plot+Sample, data=Hak.data, std.error)
TC.SE<-aggregate(Total.C..mmol.gdw~Plot+Sample, data=Hak.data, std.error)
CN.SE<-aggregate(C.N~Plot+Sample, data=Hak.data, std.error)
d13C.n<-aggregate(d13C~Plot+Sample, data=Hak.data, length)
d15N.n<-aggregate(d15N~Plot+Sample, data=Hak.data, length)
mean.data<- cbind(d13C.m, d15N.m[3], TN.m[3], TC.m[3], d13C.SE[3], d15N.SE[3], TN.SE[3], TC.SE[3], CN.m[3], CN.SE[3])
colnames(mean.data)<-c("Plot", "Sample", "d13C.mean", "d15N.mean", "TN.mean", "TC.mean", "d15N.SE", "d13C.SE", "TN.SE", "TC.SE", "CN.mean", "CN.SE")
# rename Rubus species for plotting
levels(mean.data$Sample)[levels(mean.data$Sample)=="RUBARG"] <- "Rubus spp"
levels(mean.data$Sample)[levels(mean.data$Sample)=="RUBHAW"] <- "Rubus spp"
mean.data$Sample<-factor(mean.data$Sample, levels=c("Soil", "Koa", "Rubus spp"))
### make mean plots
pd <- position_dodge(0.7) #offset for error bars
### d13C
d13C.plot<-ggplot(data=mean.data, aes(x=Sample, y=d13C.mean, fill=Plot, group=Plot))+
geom_bar(stat="identity", position = pd, width=0.7)+
geom_errorbar(aes(ymin=d13C.mean-d13C.SE, ymax=d13C.mean+d13C.SE), size=.3, width=0, position=pd) +
scale_fill_manual(values=c('gray65', "gray89")) +
xlab("Sample Type") +
annotate(geom="text", x=1, y=-27.5, label=expression(paste(bold("*"))), size=4) +
annotate(geom="text", x=2, y=-33, label=expression(paste(bold("*"))), size=4) +
annotate(geom="text", x=3, y=-33, label=expression(paste(bold("*"))), size=4) +
theme(text = element_text(size=8)) +
ylab(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))) + BW.back2
### d15N
d15N.plot<-ggplot(data=mean.data, aes(x=Sample, y=d15N.mean, fill=Plot, group=Plot))+
geom_bar(stat="identity", position = pd, width=0.7)+ ylim(c(0,8)) +
geom_errorbar(aes(ymin=d15N.mean-d15N.SE, ymax=d15N.mean+d15N.SE), size=.3, width=0, position=pd) +
scale_fill_manual(values=c('gray65', "gray89")) +
xlab("Sample Type") +
theme(text = element_text(size=8)) +
annotate(geom="text", x=2, y=2.8, label=expression(paste(bold("*"))), size=4) +
annotate(geom="text", x=3, y=3.3, label=expression(paste(bold("*"))), size=4) +
ylab(expression(paste(delta^{15}, N, " (\u2030, air)"))) + BW.back2
### total Carbon
TC.plot<-ggplot(data=mean.data, aes(x=Sample, y=TC.mean, fill=Plot, group=Plot))+
geom_bar(stat="identity", position = pd, width=0.7)+ ylim(c(0, 60)) +
geom_errorbar(aes(ymin=TC.mean-TC.SE, ymax=TC.mean+TC.SE), size=.3, width=0, position=pd) +
scale_fill_manual(values=c('gray65', "gray89")) +
xlab("Sample Type") +
theme(text = element_text(size=8)) +
annotate(geom="text", x=1, y=25, label=expression(paste(bold("*"))), size=4) +
annotate(geom="text", x=3, y=42, label=expression(paste(bold("*"))), size=4) +
ylab("Total Carbon (mmol/gdw)") + BW.back2
### total Nitrogen
TN.plot<-ggplot(data=mean.data, aes(x=Sample, y=TN.mean, fill=Plot, group=Plot))+
geom_bar(stat="identity", position = pd, width=0.7)+
geom_errorbar(aes(ymin=TN.mean-TN.SE, ymax=TN.mean+TN.SE), size=.3, width=0, position=pd) +
scale_fill_manual(values=c('gray65', "gray89")) +
xlab("Sample Type") +
theme(text = element_text(size=8)) +
ylab("Total Nitrogen (mmol/gdw)") + BW.back2
# get the legend, # create some space to the left of the legend
bar.legend <- get_legend(
TN.plot + theme(legend.box.margin = margin(0, 0, 0, 12)) +
theme(legend.key.size = unit(0.3, "cm")))
# figures together
bar.plots<-plot_grid(d15N.plot + theme(legend.position = "none"),
d13C.plot + theme(legend.position = "none"),
TN.plot + theme(legend.position = "none"),
TC.plot + theme(legend.position = "none"),
nrow=1, ncol=4, labels=c('a', 'b', 'c', 'd'),
label_size=8, hjust=-1, vjust= 3)
plot_grid(bar.plots, bar.legend, rel_widths = c(8, 1)) # legend 1/8 size as first obj.
dev.copy(pdf, "figures/Fig 2.mean.data.sample.pdf", width=7, height=3, encod="MacRoman")
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
###ISO-KRIG Give the data a “krig-over”.
First, make krigs of the full plot and all data layers
# load data
krig<-read.csv("data/krig.matrix.csv") # matrix of points for isoscape grid
# merge kriging data with isotope data
data.merge<-as.data.frame(join_all(list(Hak.data, krig), by = "Position.point", type='full'))
data.merge<-data.merge[c(-167:-194),] # drop NAs where no samples taken
scape.data<-data.merge
write.csv(scape.data, "output/scape.data.csv")
#################
####d15N isoscape
setup
################# d15N krigs
################# make site maps based on bubble plots
# just AK site
krig.AK<- scape.data %>% filter(Plot=="AK")
krig.AK$Sample<-droplevels(krig.AK$Sample)
# modify one point slightly due to overlap of coordinate system
krig.AK$y.matrix[24]=2.4; krig.AK$y.meter[24]=2.4
krig.AK$x.matrix[24]=5.4; krig.AK$x.meter[24]=5.4
AK.map<-krig.AK %>% as.data.frame %>%
ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d15N, color=Sample), alpha=3/4) +
ggtitle("AK plot--d15N (permil)") + coord_equal() + theme_bw()
# just RK site
krig.RK<- scape.data %>% filter(Plot=="RK")
krig.RK$Sample<-droplevels(krig.RK$Sample)
RK.map<-krig.RK %>% as.data.frame %>%
ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d15N, color=Sample), alpha=3/4) +
ggtitle("RK plot--d15N (permil)") + coord_equal() + theme_bw()
plot_grid(RK.map, AK.map, ncol=2, nrow=1)
dev.copy(pdf, "figures/dot.scape.pdf", width=9, height=4, encod="MacRoman")
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
#####################
#####################
AK d15N isoscape
## Krig AK sites
coordinates(krig.AK)<- ~x.meter + y.meter
# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=36, by=0.05)
Rows=seq(from=-1, to=21, by=0.05)
Grid.AK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.AK) <- ~ x+y
gridded(Grid.AK) <- TRUE # Plot the grid and points
## expand binding box
#krig.AK@bbox # current data binding box
x<-c(-1, 36)
y<-c(-1, 21)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
## min max
## x -1 36
## y -1 21
krig.AK@bbox<-bbox(S) # expanded binding box for data
Grid.AK@bbox<-bbox(S) # expand for plot corner
# autokrige
d15N.AK.auto <- autoKrige(d15N ~ 1, krig.AK, Grid.AK) # ordinary kriging
## [using ordinary kriging]
plot(d15N.AK.auto, sp.layout = list(pts = list("sp.points", krig.AK)))
#d15N.AK.auto$krige_output
# make to dataframes for lm, combine
AK.pred.N <- d15N.AK.auto[1] %>% as.data.frame() # model predictions
AK.df.N <- krig.AK %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.AK.df.N<- left_join(AK.df.N, AK.pred.N, by = c("x.meter"="krige_output.x")) # longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.AK.df.N)
summary(mod) #R squared = 0.21
##
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.AK.df.N)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6314 -0.4768 0.0794 0.4671 4.0753
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.818896 0.009430 404.99 <2e-16 ***
## d15N 0.184962 0.001937 95.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8717 on 35278 degrees of freedom
## Multiple R-squared: 0.2054, Adjusted R-squared: 0.2054
## F-statistic: 9119 on 1 and 35278 DF, p-value: < 2.2e-16
mean(pred.AK.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 4.4 d15N
## [1] 4.602741
mean(pred.AK.df.N$krige_output.var1.stdev, na.rm=T) # 2.1
## [1] 2.088436
# inspect model
# plot(pred.AK.df.N$krige_output.var1.pred ~ pred.AK.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
### map in automap
#automapPlot(d15N.AK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.AK, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# note that the prediction "var1.pred" component of the krig is a SpatialPixelsDataFrame. Need to call the krig, then the component df("krige_output"), then the output column ("var1.pred") to plot in spplot
# variance
# spplot(d15N.AK.auto$krige_output,"var1.stdev")
#predicted krige
my.palette2 <- brewer.pal(n = 6, name = "BuGn")
my.palette3 <- colorRampPalette(c("firebrick1", "darkseagreen1", "azure1"))(4)
col.scheme.N <- colorRampPalette(my.palette3)
plotd15N.krig.AK<-spplot(d15N.AK.auto$krige_output["var1.pred"], col.regions=col.scheme.N,
sp.layout = list("sp.points", krig.AK, pch=c(16,17,1)[krig.AK$Sample],
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{15}, N, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)
###########################
RK d15N isoscape
###########################
# RK Site Krig
###########################
## Krig RK site
krig.RK<- scape.data %>% filter(Plot=="RK")
krig.RK$Sample<-droplevels(krig.RK$Sample)
coordinates(krig.RK)<- ~x.meter + y.meter
Columns=seq(from=-1, to=21, by=0.05)
Rows=seq(from=-1, to=36, by=0.05)
Grid.RK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.RK) <- ~ x+y
gridded(Grid.RK) <- TRUE # Plot the grid and points
## expand binding box
#krig.RK@bbox # current data binding box
x<-c(-1, 21)
y<-c(-1, 36)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
## min max
## x -1 21
## y -1 36
Grid.RK@bbox<-bbox(S) # expand for plot corner
# autokrige
# 'fix.values' = The items describe the fixed value for the nugget (y-intercept), range (of interprolation) and sill(asymptote) respectively. Setting the value to NA means that the value is not fixed. Is passed on to autofitVariogram.
# The nugget is the y-intercept of the variogram. In practical terms, the nugget represents the small-scale variability of the data. A portion of that short range variability can be the result of measurement error.
# The range is the distance after which the variogram levels off. The physical meaning of the range is that pairs of points that are this distance or greater apart are not spatially correlated.
# The sill is the total variance contribution, or the maximum variability between pairs of points.
# increasing the range gave better fit than auto
d15N.RK.auto <- autoKrige(d15N ~ 1, krig.RK, Grid.RK, fix.value=c(NA, 2, NA)) # ordinary kriging
## [using ordinary kriging]
plot(d15N.RK.auto, sp.layout = list(pts = list("sp.points", krig.RK)))
#d15N.RK.auto$krige_output
# make to dataframes for lm, combine
RK.pred.N <- d15N.RK.auto[1] %>% as.data.frame() # model predictions
RK.df.N <- krig.RK %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
# longest axis here is y
pred.RK.df.N<- left_join(RK.df.N, RK.pred.N, by = c("x.meter"="krige_output.x"))
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.RK.df.N)
summary(mod) #R squared = 0.21
##
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.RK.df.N)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7946 -0.4047 0.0087 0.3657 3.4740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.026386 0.007487 537.8 <2e-16 ***
## d15N 0.184935 0.001434 128.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7525 on 62983 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2089, Adjusted R-squared: 0.2088
## F-statistic: 1.663e+04 on 1 and 62983 DF, p-value: < 2.2e-16
mean(pred.RK.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 5.0 d15N
## [1] 4.911052
mean(pred.RK.df.N$krige_output.var1.stdev, na.rm=T) # 1.9
## [1] 1.90181
# inspect model
# plot(pred.RK.df.N$krige_output.var1.pred ~ pred.RK.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
### map in automap
#automapPlot(d15N.RK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.RK), col.regions=my.palette)
# variance
# spplot(d15N.RK.auto$krige_output,"var1.stdev")
#predicted krige
plotd15N.krig.RK<-spplot(d15N.RK.auto$krige_output["var1.pred"], col.regions=col.scheme.N,
sp.layout = list("sp.points", krig.RK, pch=c(16,17,1)[krig.RK$Sample],
col="black", cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{15}, N, " - RK (Remnant Forest)")), cex=0.7))
#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
print(x)
}
### combined Ak-Rk krig d15N plot
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd15N.krig.AK, plotd15N.krig.RK, ncol=2, nrow=1, widths = c(1.2,1))
dev.copy(png, "figures/RKAK.d15N.krig3.png", width = 8, height = 4, units = 'in', res = 300)
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
setup
################# d13C krigs
################# make site maps based on bubble plots
AK.map<-krig.AK %>% as.data.frame %>%
ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d13C, color=Sample), alpha=3/4) +
ggtitle("AK plot--d13C (permil)") + coord_equal() + theme_bw()
RK.map<-krig.RK %>% as.data.frame %>%
ggplot(aes(x.meter, y.meter, group=Sample)) + geom_point(aes(size=d13C, color=Sample), alpha=3/4) +
ggtitle("RK plot--d13C (permil)") + coord_equal() + theme_bw()
plot_grid(RK.map, AK.map, ncol=2, nrow=1)
dev.copy(pdf, "figures/dot.scape.d13C.pdf", width=9, height=4, encod="MacRoman")
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
#####################
#####################
AK d13C isoscape
# autokrige
d13C.AK.auto <- autoKrige(d13C ~ 1, krig.AK, Grid.AK) # ordinary kriging
## [using ordinary kriging]
plot(d13C.AK.auto, sp.layout = list(pts = list("sp.points", krig.AK)))
#d13C.AK.auto$krige_output
# make to dataframes for lm, combine
AK.pred.C <- d13C.AK.auto[1] %>% as.data.frame() # model predictions
AK.df.C <- krig.AK %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.AK.df.C<- left_join(AK.df.C, AK.pred.C, by = c("x.meter"="krige_output.x")) # longest axis is x
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.AK.df.C)
summary(mod) #R squared = 0.24
##
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.AK.df.C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3830 -0.4285 0.0828 0.5165 4.2039
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -21.5295 0.0423 -509.0 <2e-16 ***
## d13C 0.1677 0.0016 104.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7706 on 35278 degrees of freedom
## Multiple R-squared: 0.2375, Adjusted R-squared: 0.2374
## F-statistic: 1.099e+04 on 1 and 35278 DF, p-value: < 2.2e-16
mean(pred.AK.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions -25.9 d13C
## [1] -25.9423
mean(pred.AK.df.C$krige_output.var1.stdev, na.rm=T) # 2.3
## [1] 2.34059
# inspect model
# plot(pred.AK.df.C$krige_output.var1.pred ~ pred.AK.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d13C.AK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.AK, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d13C.AK.auto$krige_output,"var1.stdev")
#predicted krige
my.palette4 <- colorRampPalette(c("firebrick1", "paleturquoise2", "azure1"))(4)
col.scheme.C <- colorRampPalette(my.palette4)
plotd13C.krig.AK<-spplot(d13C.AK.auto$krige_output["var1.pred"], col.regions=col.scheme.C,
sp.layout = list("sp.points", krig.AK, pch=c(16,17,1)[krig.AK$Sample],
col="black", cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{13}, C, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)
###########################
RK d13C isoscape
###########################
# RK Site Krig
###########################
## Krig RK site
# autokrige
d13C.RK.auto <- autoKrige(d13C ~ 1, krig.RK, Grid.RK, fix.value=c(NA, 1, NA)) # ordinary kriging
## [using ordinary kriging]
plot(d13C.RK.auto, sp.layout = list(pts = list("sp.points", krig.RK)))
#d13C.RK.auto$krige_output
# make to dataframes for lm, combine
RK.pred.C <- d13C.RK.auto[1] %>% as.data.frame() # model predictions
RK.df.C <- krig.RK %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
# longest axis here is y
pred.RK.df.C<- left_join(RK.df.C, RK.pred.C, by = c("x.meter"="krige_output.x"))
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.RK.df.C)
summary(mod) #R squared = 0.33
##
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.RK.df.C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7059 -0.7070 0.0331 0.7676 4.2973
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.013291 0.043240 -462.8 <2e-16 ***
## d13C 0.271933 0.001551 175.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.093 on 62983 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.3279, Adjusted R-squared: 0.3279
## F-statistic: 3.073e+04 on 1 and 62983 DF, p-value: < 2.2e-16
mean(pred.RK.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions -27.6 d13C
## [1] -27.55484
mean(pred.RK.df.C$krige_output.var1.stdev, na.rm=T) # 2.8
## [1] 2.50984
# inspect model
# plot(pred.RK.df.C$krige_output.var1.pred ~ pred.RK.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# map in automap
#automapPlot(d13C.RK.auto$krige_output, "var1.pred", sp.layout = list("sp.points", krig.RK), col.regions=my.palette)
# variance
# spplot(d13C.RK.auto$krige_output,"var1.stdev")
#predicted krige
plotd13C.krig.RK<-spplot(d13C.RK.auto$krige_output["var1.pred"], col.regions=col.scheme.C,
sp.layout = list("sp.points", krig.RK, pch=c(16,17,1)[krig.RK$Sample],
col="black", cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{13}, C, " - RK (Remnant Forest)")), cex=0.7))
#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
print(x)
}
### combined Ak-Rk krig d13C plot
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd13C.krig.AK,plotd13C.krig.RK, ncol=2, nrow=1, widths = c(1.2,1))
dev.copy(png, "figures/RKAK.d13C.krig2.png", width = 8, height = 4, units = 'in', res = 300)
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
Density plots using predictions of N and C values for all sample types
###################
# make density plots using predictions
# combine the predictions for each plot
AK.df.out.N<-as.data.frame(AK.pred.N$krige_output.var1.pred); AK.df.out.N$Plot<-"AK"
RK.df.out.N<-as.data.frame(RK.pred.N$krige_output.var1.pred); RK.df.out.N$Plot<-"RK";
colnames(AK.df.out.N)<-c("d15N.pred", "Plot")
colnames(RK.df.out.N)<-c("d15N.pred", "Plot")
# new dataframe
pred.dat.N<-rbind(AK.df.out.N, RK.df.out.N)
pred.dat.N$Plot<-as.factor(pred.dat.N$Plot)
pred.dat.N %>%
group_by(Plot) %>%
summarise(no_rows = length(Plot))
## no_rows
## 1 653562
#man whitney for significance (not normal data)
mwu(pred.dat.N, d15N.pred, Plot) # signif difference in predictions (p<0.001)
##
## # Mann-Whitney-U-Test
##
## Groups 1 = AK (n = 326781) | 2 = RK (n = 326781):
## U = 83294984480.000, W = 29901910109.000, p < 0.001, Z = -308.086
## effect-size r = 0.381
## rank-mean(1) = 254895.43
## rank-mean(2) = 398667.57
# Density plot
#hist(pred.dat.N$d15N.pred[(pred.dat.N$Plot=="AK")], ylim=c(0,150000), xlab="d15N-predict")
#hist(pred.dat.N$d15N.pred[(pred.dat.N$Plot=="RK")], ylim=c(0,150000), xlab="d15N-predict")
# plot.means
predict.mean.N <- ddply(pred.dat.N, "Plot", summarise, grp.mean=mean(d15N.pred, na.rm=TRUE))
density.predict.N<-ggplot(pred.dat.N, aes(x=d15N.pred)) +
geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+
xlab(expression(paste(delta^{15}, N["predicted"], sp="")))+
scale_color_manual(values=c("#E69F00", "#56B4E9")) +
scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
geom_vline(data=predict.mean.N, aes(xintercept=grp.mean, color=Plot),
linetype="dashed", lwd=0.2) + BW.back2
######## Carbon
###################
# make density plots using predictions
# combine the predictions for each plot
AK.df.out.C<-as.data.frame(AK.pred.C$krige_output.var1.pred); AK.df.out.C$Plot<-"AK"
RK.df.out.C<-as.data.frame(RK.pred.C$krige_output.var1.pred); RK.df.out.C$Plot<-"RK";
colnames(AK.df.out.C)<-c("d13C.pred", "Plot")
colnames(RK.df.out.C)<-c("d13C.pred", "Plot")
# new dataframe
pred.dat.C<-rbind(AK.df.out.C, RK.df.out.C)
pred.dat.C$Plot<-as.factor(pred.dat.C$Plot)
pred.dat.C %>%
group_by(Plot) %>%
summarise(no_rows = length(Plot))
## no_rows
## 1 653562
#man whitney for significance (not normal data)
mwu(pred.dat.C, d13C.pred, Plot) # signif difference in predictions (p<0.001)
##
## # Mann-Whitney-U-Test
##
## Groups 1 = AK (n = 326781) | 2 = RK (n = 326781):
## U = 147098098591.500, W = 93705024220.500, p < 0.001, Z = 528.598
## effect-size r = 0.654
## rank-mean(1) = 450142.75
## rank-mean(2) = 203420.25
# Density plot
#hist(pred.dat.C$d13C.pred[(pred.dat.C$Plot=="AK")], ylim=c(0,150000), xlab="d13C-predict")
#hist(pred.dat.C$d13C.pred[(pred.dat.C$Plot=="RK")], ylim=c(0,150000), xlab="d13C-predict")
# plot.means
predict.mean.C <- ddply(pred.dat.C, "Plot", summarise, grp.mean=mean(d13C.pred, na.rm=TRUE))
density.predict.C<-ggplot(pred.dat.C, aes(x=d13C.pred)) +
geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+
xlab(expression(paste(delta^{13}, C["predicted"], sp="")))+
scale_color_manual(values=c("#E69F00", "#56B4E9")) +
scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
geom_vline(data=predict.mean.C, aes(xintercept=grp.mean, color=Plot),
linetype="dashed", lwd=0.2) + BW.back2
# get the legend, # create some space to the left of the legend
density.legend <- get_legend(
density.predict.C + theme(legend.box.margin = margin(0, 0, 0, 12)) +
theme(legend.key.size = unit(0.3, "cm")))
gridExtra::grid.arrange(density.predict.N + theme(legend.position = "none"),
density.predict.C + theme(legend.position = "none"),
density.legend,
ncol=3, nrow=1, widths = c(1,1,0.4))
dev.copy(pdf, "figures/Fig S2ab.mod.pred.NC.pdf", width = 6, height = 2.5)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
Subset data and only examine SOIL.
- Restored forest (AK) and Remnant Forest (RK)
###### AK site d15N
## d15N SOIL plot
krig.AK<- scape.data %>% filter(Plot=="AK") # just AK
AK.soil<- krig.AK[(krig.AK$Sample=="Soil"),] #just soil
AK.soil$Sample<-droplevels(AK.soil$Sample)
## Krig AK sites
coordinates(AK.soil)<- ~x.meter + y.meter
# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=36, by=0.05)
Rows=seq(from=-1, to=21, by=0.05)
Grid.AK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.AK) <- ~ x+y
gridded(Grid.AK) <- TRUE # Plot the grid and points
## expand binding box
#AK.soil@bbox # current data binding box
x<-c(-1, 36)
y<-c(-1, 21)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
## min max
## x -1 36
## y -1 21
AK.soil@bbox<-bbox(S) # expanded binding box for data
Grid.AK@bbox<-bbox(S) # expand for plot corner
# autokrige
d15N.AK.soil <- autoKrige(d15N ~ 1, AK.soil, Grid.AK) # ordinary kriging
## [using ordinary kriging]
plot(d15N.AK.soil, sp.layout = list(pts = list("sp.points", AK.soil)))
#d15N.AK.soil$krige_output
# make to dataframes for lm, combine
AK.soil.pred.N <- d15N.AK.soil[1] %>% as.data.frame() # model predictions
AK.soil.df.N <- AK.soil %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.AK.soil.df.N<- left_join(AK.soil.df.N, AK.soil.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.AK.soil.df.N)
summary(mod) #R squared = 0.015
##
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.AK.soil.df.N)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3313 -0.5382 -0.0558 0.6719 3.0052
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.31500 0.03730 142.49 <2e-16 ***
## d15N 0.11105 0.00609 18.24 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9752 on 21166 degrees of freedom
## Multiple R-squared: 0.01547, Adjusted R-squared: 0.01542
## F-statistic: 332.6 on 1 and 21166 DF, p-value: < 2.2e-16
mean(pred.AK.soil.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 6.0 d15N
## [1] 5.984172
mean(pred.AK.soil.df.N$krige_output.var1.stdev, na.rm=T) # 0.35
## [1] 0.3514554
# inspect model
# plot(pred.AK.soil.df.N$krige_output.var1.pred ~ pred.AK.soil.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d15N.AK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", AK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d15N.AK.soil$krige_output,"var1.stdev")
#predicted krige
plotd15N.krig.AK.soil<-spplot(d15N.AK.soil$krige_output["var1.pred"], col.regions=col.scheme.N,
sp.layout = list("sp.points", AK.soil, pch=1,
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{15}, N, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)
###########################
# RK Site Krig
###########################
## d15N SOIL plot
krig.RK<- scape.data %>% filter(Plot=="RK") # just AK
RK.soil<- krig.RK[(krig.RK$Sample=="Soil"),] #just soil
RK.soil$Sample<-droplevels(RK.soil$Sample)
## Krig RK sites
coordinates(RK.soil)<- ~x.meter + y.meter
# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=21, by=0.05)
Rows=seq(from=-1, to=36, by=0.05)
Grid.RK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.RK) <- ~ x+y
gridded(Grid.RK) <- TRUE # Plot the grid and points
## expand binding box
#RK.soil@bbox # current data binding box
x<-c(-1, 21)
y<-c(-1, 36)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
## min max
## x -1 21
## y -1 36
RK.soil@bbox<-bbox(S) # expanded binding box for data
Grid.RK@bbox<-bbox(S) # expand for plot corner
# autokrige
d15N.RK.soil <- autoKrige(d15N ~ 1, RK.soil, Grid.RK) # ordinary kriging
## [using ordinary kriging]
plot(d15N.RK.soil, sp.layout = list(pts = list("sp.points", RK.soil)))
#d15N.RK.soil$krige_output
# make to dataframes for lm, combine
RK.soil.pred.N <- d15N.RK.soil[1] %>% as.data.frame() # model predictions
RK.soil.df.N <- RK.soil %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.RK.soil.df.N<- left_join(RK.soil.df.N, RK.soil.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.RK.soil.df.N)
summary(mod) #R squared = 0.01
##
## Call:
## lm(formula = krige_output.var1.pred ~ d15N, data = pred.RK.soil.df.N)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.87178 -0.34756 0.02653 0.50770 1.97596
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.124328 0.022592 271.08 <2e-16 ***
## d15N 0.052800 0.003492 15.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.708 on 35566 degrees of freedom
## Multiple R-squared: 0.006386, Adjusted R-squared: 0.006358
## F-statistic: 228.6 on 1 and 35566 DF, p-value: < 2.2e-16
mean(pred.RK.soil.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 6.4 d15N
## [1] 6.461148
mean(pred.RK.soil.df.N$krige_output.var1.stdev, na.rm=T) # 0.8
## [1] 0.7599088
# inspect model
# plot(pred.RK.soil.df.N$krige_output.var1.pred ~ pred.RK.soil.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d15N.RK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", RK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d15N.RK.soil$krige_output,"var1.stdev")
#predicted krige
plotd15N.krig.RK.soil<-spplot(d15N.RK.soil$krige_output["var1.pred"], col.regions=col.scheme.N,
sp.layout = list("sp.points", RK.soil, pch=1,
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{15}, N, " - RK (Remnant Forest)")), cex=0.8), colorkey=TRUE)
#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
print(x)
}
### combined Ak-Rk krig d15N plot
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd15N.krig.AK.soil, plotd15N.krig.RK.soil, ncol=2, nrow=1, widths = c(1.2,1))
dev.copy(png, "figures/RKAK.d15N.soils.png", width = 8, height = 4, units = 'in', res = 300)
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
###### AK site d13C
# autokrige
d13C.AK.soil <- autoKrige(d13C ~ 1, AK.soil, Grid.AK) # ordinary kriging
## [using ordinary kriging]
plot(d13C.AK.soil, sp.layout = list(pts = list("sp.points", AK.soil)))
#d13C.AK.soil$krige_output
# make to dataframes for lm, combine
AK.soil.pred.C <- d13C.AK.soil[1] %>% as.data.frame() # model predictions
AK.soil.df.C <- AK.soil %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.AK.soil.df.C<- left_join(AK.soil.df.C, AK.soil.pred.C, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.AK.soil.df.C)
summary(mod) #R squared = 0.008
##
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.AK.soil.df.C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.61780 -0.52385 0.02572 0.48253 2.96360
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22.734708 0.130976 -173.58 <2e-16 ***
## d13C 0.069928 0.005368 13.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8341 on 21166 degrees of freedom
## Multiple R-squared: 0.007955, Adjusted R-squared: 0.007908
## F-statistic: 169.7 on 1 and 21166 DF, p-value: < 2.2e-16
mean(pred.AK.soil.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions -24.4 d13C
## [1] -24.43935
mean(pred.AK.soil.df.C$krige_output.var1.stdev, na.rm=T) # 0.6
## [1] 0.6411003
# inspect model
# plot(pred.AK.soil.df.C$krige_output.var1.pred ~ pred.AK.soil.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d13C.AK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", AK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d13C.AK.soil$krige_output,"var1.stdev")
#predicted krige
plotd13C.krig.AK.soil<-spplot(d13C.AK.soil$krige_output["var1.pred"], col.regions=col.scheme.C,
sp.layout = list("sp.points", AK.soil, pch=1,
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{13}, C, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)
###########################
# RK Site Krig
###########################
# autokrige
d13C.RK.soil <- autoKrige(d13C ~ 1, RK.soil, Grid.RK) # ordinary kriging
## [using ordinary kriging]
plot(d13C.RK.soil, sp.layout = list(pts = list("sp.points", RK.soil)))
#d13C.RK.soil$krige_output
# make to dataframes for lm, combine
RK.soil.pred.C <- d13C.RK.soil[1] %>% as.data.frame() # model predictions
RK.soil.df.C <- RK.soil %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.RK.soil.df.C<- left_join(RK.soil.df.C, RK.soil.pred.C, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d13C, data=pred.RK.soil.df.C)
summary(mod) #R squared = 0.16
##
## Call:
## lm(formula = krige_output.var1.pred ~ d13C, data = pred.RK.soil.df.C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.87389 -0.51641 -0.01577 0.55010 2.42974
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22.840353 0.093601 -244.02 <2e-16 ***
## d13C 0.097864 0.003685 26.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6794 on 35566 degrees of freedom
## Multiple R-squared: 0.01944, Adjusted R-squared: 0.01941
## F-statistic: 705.2 on 1 and 35566 DF, p-value: < 2.2e-16
mean(pred.RK.soil.df.C$krige_output.var1.pred, na.rm=T) # mean of predictions 4.4 d13C
## [1] -25.32413
mean(pred.RK.soil.df.C$krige_output.var1.stdev, na.rm=T) # 2.1
## [1] 0.6704758
# inspect model
# plot(pred.RK.soil.df.C$krige_output.var1.pred ~ pred.RK.soil.df.C$d13C, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d13C.RK.soil$krige_output, "var1.pred", sp.layout = list("sp.points", RK.soil, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d13C.RK.soil$krige_output,"var1.stdev")
#predicted krige
plotd13C.krig.RK.soil<-spplot(d13C.RK.soil$krige_output["var1.pred"], col.regions=col.scheme.C,
sp.layout = list("sp.points", RK.soil, pch=1,
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{13}, C, " - RK (Remnant Forest)")), cex=0.8), colorkey=TRUE)
#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
print(x)
}
### combined Ak-Rk krig d13C plot
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd13C.krig.AK.soil, plotd13C.krig.RK.soil, ncol=2, nrow=1, widths = c(1.2,1))
dev.copy(png, "figures/RKAK.d13C.soils.png", width = 8, height = 4, units = 'in', res = 300)
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
Density plots for soil isoscapes
# make density plots using predictions
# combine the predictions for each plot
############## Nitrogen plot
AK.df.soil.N<-as.data.frame(pred.AK.soil.df.N$krige_output.var1.pred); AK.df.soil.N$Plot<-"AK"
RK.df.soil.N<-as.data.frame(pred.RK.soil.df.N$krige_output.var1.pred); RK.df.soil.N$Plot<-"RK"
colnames(AK.df.soil.N)<-c("d15N.pred", "Plot")
colnames(RK.df.soil.N)<-c("d15N.pred", "Plot")
# new dataframe
pred.soil.N<-rbind(AK.df.soil.N, RK.df.soil.N)
pred.soil.N$Plot<-as.factor(pred.soil.N$Plot)
pred.soil.N %>%
group_by(Plot) %>%
summarise(no_rows = length(Plot))
## no_rows
## 1 56736
#man whitney for significance (not normal data)
mwu(pred.soil.N, d15N.pred, Plot) # signif difference in predictions (p<0.001)
##
## # Mann-Whitney-U-Test
##
## Groups 1 = AK (n = 21168) | 2 = RK (n = 35568):
## U = 474341208.000, W = 250288512.000, p < 0.001, Z = -66.868
## effect-size r = 0.281
## rank-mean(1) = 22408.41
## rank-mean(2) = 31915.60
# Density plot
#hist(pred.soil.N$d15N.pred[(pred.soil.N$Plot=="AK")], ylim=c(0,10000), xlab="d15N-predict")
#hist(pred.soil.N$d15N.pred[(pred.soil.N$Plot=="RK")], ylim=c(0,10000), xlab="d15N-predict")
# plot.means
predict.mean.N <- ddply(pred.soil.N, "Plot", summarise, grp.mean=mean(d15N.pred, na.rm=TRUE))
density.predict.Nsoil<-ggplot(pred.soil.N, aes(x=d15N.pred)) +
geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+ ylim(0,1.2) +
xlab(expression(paste(delta^{15}, N["predicted"], sp="")))+
scale_color_manual(values=c("#E69F00", "#56B4E9")) +
scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
geom_vline(data=predict.mean.N, aes(xintercept=grp.mean, color=Plot),
linetype="dashed", lwd=0.2) + BW.back2
############## Carbon
AK.df.soil.C<-as.data.frame(pred.AK.soil.df.C$krige_output.var1.pred); AK.df.soil.C$Plot<-"AK"
RK.df.soil.C<-as.data.frame(pred.RK.soil.df.C$krige_output.var1.pred); RK.df.soil.C$Plot<-"RK"
colnames(AK.df.soil.C)<-c("d13C.pred", "Plot")
colnames(RK.df.soil.C)<-c("d13C.pred", "Plot")
# new dataframe
pred.soil.C<-rbind(AK.df.soil.C, RK.df.soil.C)
pred.soil.C$Plot<-as.factor(pred.soil.C$Plot)
pred.soil.C %>%
group_by(Plot) %>%
summarise(no_rows = length(Plot))
## no_rows
## 1 56736
#man whitney for significance (not normal data)
mwu(pred.soil.C, d13C.pred, Plot) # signif difference in predictions (p<0.001)
##
## # Mann-Whitney-U-Test
##
## Groups 1 = AK (n = 21168) | 2 = RK (n = 35568):
## U = 819376896.000, W = 595324200.000, p < 0.001, Z = 116.006
## effect-size r = 0.487
## rank-mean(1) = 38708.28
## rank-mean(2) = 22214.87
# Density plot
#hist(pred.soil.C$d13C.pred[(pred.soil.C$Plot=="AK")], ylim=c(0,10000), xlab="d13C-predict")
#hist(pred.soil.C$d13C.pred[(pred.soil.C$Plot=="RK")], ylim=c(0,10000), xlab="d13C-predict")
# plot.means
predict.mean.C <- ddply(pred.soil.C, "Plot", summarise, grp.mean=mean(d13C.pred, na.rm=TRUE))
density.predict.Csoil<-ggplot(pred.soil.C, aes(x=d13C.pred)) +
geom_density(aes(color=Plot, fill=Plot), alpha=0.5)+ ylim(0,0.8) +
xlab(expression(paste(delta^{13}, C["predicted"], sp="")))+
scale_color_manual(values=c("#E69F00", "#56B4E9")) +
scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
geom_vline(data=predict.mean.C, aes(xintercept=grp.mean, color=Plot),
linetype="dashed", lwd=0.2) + BW.back2
# get the legend, # create some space to the left of the legend
density.legend <- get_legend(
density.predict.Nsoil + theme(legend.box.margin = margin(0, 0, 0, 12)) +
theme(legend.key.size = unit(0.3, "cm")))
gridExtra::grid.arrange(density.predict.Nsoil + theme(legend.position = "none"),
density.predict.Csoil + theme(legend.position = "none"),
density.legend,
ncol=3, nrow=1, widths = c(1,1,0.3))
dev.copy(pdf, "figures/Fig Sx.mod.pred.pdf", width = 5, height = 2.5)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
#### d15N Plants isoscape
Subset data and only examine PLANTS
- *Restored forest (AK)* and *Remnant Forest (RK)*
###### AK site d15N
## d15N plants plot
krig.AK<- scape.data %>% filter(Plot=="AK") # just AK
AK.plants<- krig.AK[!(krig.AK$Sample=="Soil"),] #just plants
AK.plants$Sample<-droplevels(AK.plants$Sample)
## Krig AK sites
coordinates(AK.plants)<- ~x.meter + y.meter
# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=36, by=0.05)
Rows=seq(from=-1, to=21, by=0.05)
Grid.AK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.AK) <- ~ x+y
gridded(Grid.AK) <- TRUE # Plot the grid and points
## expand binding box
#AK.plants@bbox # current data binding box
x<-c(-1, 36)
y<-c(-1, 21)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
AK.plants@bbox<-bbox(S) # expanded binding box for data
Grid.AK@bbox<-bbox(S) # expand for plot corner
# autokrige
d15N.AK.plants <- autoKrige(d15N ~ 1, AK.plants, Grid.AK, fix.value=c(0.2, 3, 0.6)) # ordinary kriging
plot(d15N.AK.plants, sp.layout = list(pts = list("sp.points", AK.plants)))
#d15N.AK.plants$krige_output
# make to dataframes for lm, combine
AK.plants.pred.N <- d15N.AK.plants[1] %>% as.data.frame() # model predictions
AK.plants.df.N <- AK.plants %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.AK.plants.df.N<- left_join(AK.plants.df.N, AK.plants.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.AK.plants.df.N)
summary(mod) #R squared = 0.19
mean(pred.AK.plants.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 1.6 d15N
mean(pred.AK.plants.df.N$krige_output.var1.stdev, na.rm=T) # 0.7
# inspect model
# plot(pred.AK.plants.df.N$krige_output.var1.pred ~ pred.AK.plants.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d15N.AK.plants$krige_output, "var1.pred", sp.layout = list("sp.points", AK.plants, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d15N.AK.plants$krige_output,"var1.stdev")
#predicted krige
plotd15N.krig.AK.plants<-spplot(d15N.AK.plants$krige_output["var1.pred"], col.regions=col.scheme.N,
sp.layout = list("sp.points", AK.plants, pch=c(16,17)[AK.plants$Sample],
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{15}, N, " - AK (Restored Forest)")), cex=0.8), colorkey=FALSE)
plotd15N.krig.AK.plants
###########################
# RK Site Krig
###########################
## d15N plants plot
krig.RK<- scape.data %>% filter(Plot=="RK") # just AK
RK.plants<- krig.RK[!(krig.RK$Sample=="Soil"),] #just plants
RK.plants$Sample<-droplevels(RK.plants$Sample)
## Krig RK sites
coordinates(RK.plants)<- ~x.meter + y.meter
# Create a grid of "Pixels" using x as columns and y as rows
# And the rows and columns of pixels:
Columns=seq(from=-1, to=21, by=0.05)
Rows=seq(from=-1, to=36, by=0.05)
Grid.RK <- expand.grid(x=Columns,y=Rows)
coordinates(Grid.RK) <- ~ x+y
gridded(Grid.RK) <- TRUE # Plot the grid and points
## expand binding box
#RK.plants@bbox # current data binding box
x<-c(-1, 21)
y<-c(-1, 36)
xy<-cbind(x,y)
S<-SpatialPoints(xy)
bbox(S)
RK.plants@bbox<-bbox(S) # expanded binding box for data
Grid.RK@bbox<-bbox(S) # expand for plot corner
# autokrige
d15N.RK.plants <- autoKrige(d15N ~ 1, RK.plants, Grid.RK) # ordinary kriging
plot(d15N.RK.plants, sp.layout = list(pts = list("sp.points", RK.plants)))
#d15N.RK.plants$krige_output
# make to dataframes for lm, combine
RK.plants.pred.N <- d15N.RK.plants[1] %>% as.data.frame() # model predictions
RK.plants.df.N <- RK.plants %>% as.data.frame() # dataframe
# make dataframe with data and predictions (krig) to assess fit
pred.RK.plants.df.N<- left_join(RK.plants.df.N, RK.plants.pred.N, by = c("x.meter"="krige_output.x"))#longest axis is x
mod <- lm(krige_output.var1.pred ~ d15N, data=pred.RK.plants.df.N)
summary(mod) #R squared = 0.15
mean(pred.RK.plants.df.N$krige_output.var1.pred, na.rm=T) # mean of predictions 2.8 d15N
mean(pred.RK.plants.df.N$krige_output.var1.stdev, na.rm=T) # 0.5
# inspect model
# plot(pred.RK.plants.df.N$krige_output.var1.pred ~ pred.RK.plants.df.N$d15N, xlab = "Observed", ylab = "Predicted", cex=0.2); abline(0,1, lty=2); abline(mod, col = "blue")
# Adding the sp.layout parameter shows the locations of the measurements
my.palette <- brewer.pal(n = 9, name = "RdBu")
# map in automap
#automapPlot(d15N.RK.plants$krige_output, "var1.pred", sp.layout = list("sp.points", RK.plants, pch=16, col="gray20", cex=0.5), col.regions=my.palette)
# map in ssplot
# variance
# spplot(d15N.RK.plants$krige_output,"var1.stdev")
#predicted krige
plotd15N.krig.RK.plants<-spplot(d15N.RK.plants$krige_output["var1.pred"], col.regions=col.scheme.N,
sp.layout = list("sp.points", RK.plants, pch=c(16,17)[RK.plants$Sample],
col=c("black"), cex=0.5, alpha=0.5),
main=list(label=expression(paste(delta^{15}, N, " - RK (Remnant Forest)")), cex=0.8), colorkey=TRUE)
plotd15N.krig.RK.plants
#####################
grid.draw.ggmatrix <- function(x, recording = TRUE) {
print(x)
}
### combined Ak-Rk krig d15N plot
par(mfrow=c(1,2))
gridExtra::grid.arrange(plotd15N.krig.AK.plants, plotd15N.krig.RK.plants, ncol=2, nrow=1, widths = c(1.2,1))
dev.copy(png, "figures/RKAK.d15N.plants.png", width = 8, height = 4, units = 'in', res = 300)
dev.off()
#####################
# Carbon kriging, leave alone for now.
#####################
## Krig AK site d13C
coordinates(krig.AK)<- ~x.meter + y.meter
col.scheme.C <- colorRampPalette(c('coral', 'cadetblue2', 'gray5'))(5)
Grid.AK <- spsample(krig.AK, type='regular', n=1e5)
gridded(Grid.AK) <- TRUE
## d13C plot
d13C.krig.AK <- krige(d13C ~ 1, krig.AK, Grid.AK)
plotd13C.krig.AK<-spplot(d13C.krig.AK["var1.pred"], col.regions=colorRampPalette(col.scheme.C),
main=expression(paste(delta^{13}, C, "--AK (Afforested Koa Plot)")))
## Krig RK site d13C
coordinates(krig.RK)<- ~x.meter + y.meter
GridRK <- spsample(krig.RK, type='regular', n=1e5)
gridded(GridRK) <- TRUE
## d13C plot
d13C.krig.RK <- krige(d13C ~ 1, krig.RK, GridRK)
plotd13C.krig.RK<-spplot(d13C.krig.RK["var1.pred"], col.regions=colorRampPalette(col.scheme.C),
main=expression(paste(delta^{13}, C, "--RK (Remnant Koa Plot)")))